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We perform molecular dynamics simulations of model granular systems undergoing boundary- 
driven planar shear flow in two spatial dimensions with the goal of developing a more complete 
understanding of how dense particulate systems respond to applied shear. In particular, we are 
interested in determining when these systems will possess linear velocity profiles and when they will 
develop highly localized velocity profiles in response to shear. In previous work on similar systems 
we showed that nonlinear velocity profiles form when the speed of the shearing boundary exceeds 
the speed of shear waves in the material. However, we find that nonlinear velocity profiles in these 
systems are unstable at very long times. The degree of nonlinearity slowly decreases in time; the 
velocity profiles become linear when the granular temperature and density profiles are uniform across 
the system at long times. We measure the time ij required for the velocity profiles to become linear 
and find that ti increases as a power-law with the speed of the shearing boundary and increases 
rapidly as the packing fraction approaches random close packing. We also performed simulations 
in which differences in the granular temperature across the system were maintained by vertically 
vibrating one of the boundaries during shear flow. We find that nonlinear velocity profiles form and 
are stable at long times if the difference in the granular temperature across the system exceeds a 
threshold value that is comparable to the glass transition temperature in an equilibrium system at 
the same average density. Finally, the sheared and vibrated systems form stable shear bands, or 
highly localized velocity profiles, when the applied shear stress is lowered below the yield stress of 
the static part of the system. 

PACS numbers: 83.10.Rs,83.50.Ax,45.70.Mg,64.70.Pf 



I. INTRODUCTION 

The response of fluids to applied shear is an important 
and well studied problem. When a Newtonian fluid is 
slowly sheared, for example by moving the top bound- 
ary of the system at fixed velocity u in the cc-direction at 
height y = L y relative to a stationary bottom boundary 
at y — 0, a linear velocity profile v x (y) = 7?/ is estab- 
lished, where the shear rate is 7 = u/L y . In addition, 
in the small shear rate limit the shear stress is linearly 
related to the shear rate a xy = 777, where r\ is the shear 
viscosity. This relation is often employed to measure the 
shear viscosity of simple fluids. 

However, what is the response of non-Newtonian fluids 
like granular materials to an applied shear? In contrast 
to simple liquids, it is extremely difficult to predict the 
response of dense granular media to shear because they 
are inherently out of thermal equilibrium, interact via 
frictional and enduring contacts, and possess a nonzero 
yield stress. Granular materials do not flow homoge- 
neously when they are sheared, instead, shear is often 
localized into shear bands. When this occurs, most of 
the flow is confined to a narrow, locally dilated region 
near the shearing boundary while the remainder of the 
system is nearly static. Recent experimental work in- 
vestigating the response of dense granular materials to 
shear includes studies of Couette flow in 2D [lj, in 3D 
for spherical particles 0, and as a function of parti- 
cle shape 4], studies of cyclic planar shear la], studies 



of wide shear zones formed in the bulk using a modified 
Couette geometry 6] , and studies of chute flow . 

Recent theoretical studies Q have shown that kinetic 
theory can correctly predict the velocity, granular tem- 
perature, and density profiles found in experiments of 
Couette flow in the dilute regime 9]. However, kinetic 
theory and modifications included to account for diverg- 
ing viscosity are not analytically tractable, and are un- 
likely to predict accurately the properties of dense shear 
flows. Thus, molecular dynamics simulations are often 
employed to study dense granular shear flows, for exam- 
ple in Refs. 0, 0, 0, We choose a similar plan 
of attack and perform molecular dynamics simulations 
of model frictionless dense granular systems undergoing 
boundary-driven planar shear flow in 2D. We are inter- 
ested in answering several important questions: What are 
the velocity, granular temperature, and density profiles 
as a function of the velocity of the shearing boundary? In 
particular, do highly localized velocity profiles form and, 
if so, are they stable at long times? We will investigate 
these questions in simple 2D systems composed of inelas- 
tic but frictionless particles in planar shear cells, and in 
the absence of gravity. Our intent is to understand in 
detail the time and spatial dependence of the response 
in these simple systems to shear first, and then extend 
our studies to 3D systems, systems composed of frictional 
particles, and Couette shear cells. In Sec. IVII we present 
preliminary results from these future studies. 

We previously reported that nonlinear velocity profiles 
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form in repulsive athermal systems when the velocity of 
the shearing wall exceeds the speed of shear waves in 
the material [Til ]. In the present article, we study much 
longer time scales and show that nonlinear velocity pro- 
files are unstable at long times. In addition, the gran- 
ular temperature and density profiles become uniform 
throughout the system at long times. Thus, granular 
temperature and density gradients cannot be maintained 
by planar shear flow in dense systems with dissipative but 
frictionless interactions. We measure the time ti for the 
velocity profiles to become linear and find that ti scales 
as a power-law in the velocity of the shearing boundary, 
and increases rapidly as the average density of the sys- 
tem approaches random close packing <fi rcp from above, 
where cL r re 0.84 in 2D [T^. 

To maintain a granular temperature difference across 
the system, we also studied systems that were both 
sheared and vertically vibrated. We find that if the 
difference in the granular temperature across the sys- 
tem exceeds a threshold that is comparable to the glass 
transition temperature in an equilibrium system at the 
same average density, nonlinear velocity profiles are sta- 
ble at long times. Finally, we show that shear bands, 
or highly localized velocity profiles, form in the vibrated 
and sheared systems if the shear stress is tuned below 
the yield stress of the static part of the system. 



II. METHODS 

Before discussing our results further, we will first de- 
scribe our numerical model and methods. We per- 
formed molecular dynamics simulations of purely re- 
pulsive and frictionless athermal systems undergoing 
boundary-driven planar shear flow in two spatial dimen- 
sions. The systems were composed of N/2 large and 
N/2 small particles with equal mass m and diameter 
ratio 1.4 to prevent crystallization and segregation dur- 
ing shear. The starting configurations were prepared by 
choosing an average packing fraction and random initial 
positions and then allowing the system to relax to the 
nearest local potential energy minimum [l5| using the 
conjugate gradient method |l6j |. During the quench, pe- 
riodic boundary conditions were implemented in both the 
x- and y-directions. Following the quench, particles with 
y-coordinates y > L y (y < 0) were chosen to comprise 
the top (bottom) boundary. Thus, the top and bottom 
walls were rough and amorphous. 

Shear flow in the x-direction with a shear gradient in 
the y-direction and global shear rate u/L y was created 
by moving all particles in the top wall at fixed velocity 
u in the ^-direction relative to the stationary bottom 
wall. During shear flow, periodic boundary conditions 
were imposed in the x-direction. We chose an aspect 
ratio L x /L y — 1/4 with more than 50 particles along the 
shear gradient direction to reduce finite-size effects, and 
focused on systems with packing fractions in the range 
= [0.835,0.95]. 



Both bulk particles and particles comprising the 
boundary interact via the purely repulsive harmonic 
spring potential 

where e is the characteristic energy scale of the interac- 
tion, <7jj = (<7i + aj)/2 is the average diameter of particles 
i and j, rij is their separation, and Q(x) is the Heaviside 
step function. Note that the interaction potential is zero 
when > ov,- . 

In our simulations, we employ athermal or dissipative 
dynamics with no frictional or tangential forces [17j . The 
position and velocity of particles in the bulk were ob- 
tained by solving 

d 2 r 

m ~ d ~^ = F i ~ 5 « Yl ^ ~ Vj) ■ ^ij] hj , (2) 
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where FT = -£\ dV(r tj )/d 

fij^ij: the sums over j only 
include particles that overlap i, Vi is the velocity of parti- 
cle i, and b n > is the damping coefficient. We focus on 
underdamped systems in this study; specifically, for most 
simulations we use b n = 0.0375 (coefficient of restitution 
e = 0.92). In addition, we show some results that cover 
a range of e = [0.1,0.99]. The units of length, energy, 
and time are ch osen as the small particle diameter <r, e, 
and 1/lJc = cryrnje, and all quantities were normalized 
by these. 

We also performed simulations in which the system was 
both sheared and vertically vibrated. As discussed above, 
the system was sheared by constraining the top boundary 
to move in the x-direction at fixed speed u relative to the 
bottom boundary. In addition, particles in either the top 
or bottom boundary were vibrated vertically such that 
the y-coordinates of the boundary particles varied in time 
as 

Vi = Vio + A sin (ut) , (3) 

where yio is the initial position of boundary particle i, 
A is the amplitude, and to is the angular frequency of 
the vibration. We fixed A = a/5 so that the vibrations 
do not lead to large, unphysical particle overlaps and u> 
was tuned over a range of frequencies. Note that in our 
choice of units u> is normalized by the natural frequency 
lu c of the linear spring interactions. It should also be 
pointed out that the simulations with both shear and 
vibration are not performed at constant volume but at 
fixed average volume. 

We measured several physical quantities in these simu- 
lations including the local flow velocity v x , packing frac- 
tion 4>, and velocity fluctuations or granular temperature 
(Sv y ) 2 = (Vy) — (v y ) 2 in the shear gradient direction, 
as a function of the boundary velocity u and vibration 
frequency w. To study properties as a function of the ver- 
tical distance from the fixed boundaries, we divided the 
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FIG. 1: Shear strain 7 = x/L v versus time t for two shear 
stresses a xy = 1.57 x 10~ 3 (solid line) and 1.80 x 10~ 3 (dashed 
line) applied to a system at cj> = 0.90 that was unsheared at 
t = 0. 

system into rectangular bins centered at height y and av- 
eraged the quantities over the height of the bin Ay w 2 
large particle diameters. To improve the statistics and 
to study time dependence, we also performed ensemble 
averages over at least 50 different initial configurations. 
In the discussion below, we will denote ensemble averages 
using (.). When the systems are stationary but still fluc- 
tuate in time, we also perform time averages and these 
are denoted by {.) t . 

We also monitored the shear stress a xy during the 
course of our simulations. The total shear stress in the 
bulk can be calculated using the virial expression 

<Jxy = I X! Sv i,* Sv i,V + X! r V,* F V>V ' ( 4 ) 

X V \ i i>3 ) 

where Tij x is the ^-component of = r*j — fj, Fij,y 
is the y-componcnt of the total pair force Fij including 
both conservative and damping forces, and Svi lX and 8vi >y 
measure deviations in the velocity of a particle from the 
velocity averaged over the entire system. We also com- 
pared the shear stress obtained from Eq.0|with the shear 
stress \F^\/L X on the boundaries, where F^ is the total 
force in the x-direction acting on either the top or bot- 
tom boundary. On average, the bulk shear stress in Eq.0] 
and the shear stress on the boundaries were within a few 
percent of each other. 

To make contact with the recent results on sheared 
Lennard- Jones glasses in Ref. 1 we measured the yield 
stress of the static starting configurations. To do this, 
we applied a constant horizontal force F (or shear stress 
G xy = F / L x ) to the top boundary. Results from the sim- 
ulations at constant horizontal force are shown in Fig.^ 
Initially, the system flows. However, when the shear 
stress is below the yield stress, the system finds a state 
that can sustain the applied shear stress and stops flow- 
ing. If the applied shear stress is above the yield stress, 
the system will flow indefinitely. We defined the yield 



FIG. 2: Ensemble averaged velocity (v x (y)) in the shear flow 
direction as a function of height y/L v from the stationary 
boundary for a system sheared at u = 0.364 with average 
packing fraction (f> — 0.90 and b n — 0.0375 at several different 
times t > t s ~ 470, where t s is the time required for a shear 
wave to traverse the system. The curves correspond to t = 
520 (dashed line), 1000 (dotted line), 2000 (dot-dashed line), 
and 20000 (solid line). These times correspond to large shear 
strains 7 = 2.6, 5, 10, and 100. Each curve was averaged over 
an ensemble of at least 50 different starting configurations. 

shear stress cto as the minimum shear stress above which 
the shear strain continues to increase beyond 7 = 10. 

III. RESULTS 

In this section, we report the results from two sets of 
numerical simulations of purely repulsive and frictionless 
athermal systems. Section IIII Al presents results from 
simulations of boundary-driven planar shear flow. We 
show that nonlinear velocity profiles form at short time 
scales, but they slowly evolve into linear profiles at long 
times. As the velocity profiles evolve toward linear ones, 
the local granular temperature and packing fraction be- 
come uniform. Section IlII Bl presents results from simu- 
lations of boundary-driven planar shear flow in the pres- 
ence of vertical vibrations. We find that highly nonlinear 
velocity profiles can be stabilized at long times if a suf- 
ficiently large granular temperature difference is main- 
tained across the system. 

A. Time evolution of velocity profiles 

In our prev ious studies of boundary-driven planar 
shear flow |l4j , we reported that nonlinear velocity pro- 
files form when the velocity u of the shearing boundary 
exceeds u c = u s /2, where u s is speed of shear waves in 
the system. This condition was obtained by comparing 
the time t s = 2L y /u s for a shear wave to traverse the sys- 
tem to the time t u = L y /u for the system to shear unit 
strain. A rough estimate of the speed of shear waves (at 
least in the low shear rate limit) can be obtained from 
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FIG. 3: Degree of nonlinearity (A(t)) y of the velocity profile 
as a function of time t at several different boundary velocities 
in a system at cj) = 0.90 and b n — 0.0375. We show boundary 
velocities u = 0.00727, 0.0145, 0.0364, 0.0727, 0.145, 0.364, 
u — 0.727 over two orders of magnitude; u increases from left 
to right as shown by the arrow. The dotted line {A(t)) y — 
0.05 was used to estimate the time ti required for the system 
to attain an approximately linear velocity profile. The spike in 
{A(t)) y at small u and short times occurs because the velocity 
profile can fluctuate above the linear profile. 



\J G I p, where G is the static shear modulus and p is 
the mass density. A more precise way to measure u s is 
to calculate the transverse current correlation function 
Cx(k, v) [l9| as a function of frequency v and wave num- 
ber k — 2irn/L x (n — integer) and determine the slope 
of the resulting dispersion relation v{k) [2jj. One of the 
novel aspects of our previous work was that we showed 
that underdamped systems would be extremely suscepti- 
ble to nonlinear velocity profiles near random close pack- 
ing since the critical velocity u c tends to zero at <p rcp - 
We therefore predicted that any u > would give rise to 
nonlinear velocity profiles in these systems near </> rcp - 

In these prior studies, we measured the velocity, pack- 
ing fraction, and mean-square velocity fluctuation pro- 
files after shearing the system for a strain of at least 5 
and times t > t s so that the shear stress had relaxed to 
its long-time average value. This protocol for bringing 
sheared systems to steady-state is typical in both sim- 
ulations and experiments. If t s were the only relevant 
time scale, the nonlinear velocity profiles that occur in 
boundary-driven planar shear flow when u > u c would 
be stable over long times. However, in more recent stud- 
ies, we have found that these nonlinear velocity profiles 
are not stable at long times t t s and slowly evolve 
toward linear profiles. In Fig. |3 we show the slow time 
evolution of the velocity profile in a system sheared at 
u = 0.364, <j) = 0.90, and e = 0.92. Note that the strains 
beyond which the profiles become linear are extremely 
large, 7 > 25. We will also show below that the time re- 
quired for the velocity profile to become linear increases 
as the system becomes more elastic. Most of our previ- 
ous work focused on nearly elastic systems with e « 0.98 
and timescales near t s , which made it difficult to detect 



4 


2 r 

4 


7 3 


8 - 


° 3 


6 


^ 3 


4 - 


3 


2 

-A 



-3.5 -3 -2.5 
log 10 (U-U,)/L 



-2 



FIG. 4: Time ti required for the velocity profiles to reach an 
approximately linear profile relative to the time t a required 
for a shear wave to traverse the system, as a function of the 
velocity of the shearing boundary u — Ui. u% is the velocity 
of the shearing boundary below which ti = t s . The solid 
line has slope 0.5. The system parameters are <j> = 0.90 and 
b n = 0.0375. 

the profile's slow evolution. It would be interesting to 
know whether this slow evolution can be seen in exper- 
iments on sheared granular media, or other particulate 
systems. However, few experiments have studied such 
large strain and time scales. We note that recent exper- 
iments on granular media undergoing planar cyclic 
shear do show slow evolution; however, further experi- 
ments are required. Similarly, simulations of frictional 
granular media undergoing Couette flow in 3D have re- 
ported significant time evolution of the measured velocity 
profiles |21| . 

To quantify the shape of the velocity profiles as a func- 
tion of time, we define the degree of nonlinearity as 



A(t) 



1 



{vx(y))t 



(5) 



where (v x (y))t is the time and ensemble average of the ve- 
locity in the flow direction after the evolution of the pro- 
file has ended. To measure the degree of nonlinearity, we 
calculated (A(t)) y averaged over the central part of the 
system excluding layers that are immediately adjacent to 
the top and bottom walls. (A(t)) y « corresponds to a 
nearly linear velocity profile, while (A(t)) y « 1 is highly 
nonlinear. 

In Fig. |3 we show (A(t)) y for a system with <f> — 0.90 
and e = 0.92 at several different velocities of the shearing 
boundary. At each u, the velocity profiles slowly evolve 
toward linear profiles and (A(t)) y decays to zero at long 
times. To characterize the long-time behavior, we define 
a timescale U as the time required for (A(t)) y to decay to 
0.05, which is slightly above the noise level. By definition, 
the velocity profiles are steady and linear for times t > ti. 
At large u, the velocity profiles approach the linear profile 
from below as shown in Fig. [5] At small u, {A(t)) y decays 
to zero quickly, but the profiles jump above and below the 
linear profile as t — * ti . The fact that the velocity profile 
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FIG. 5: Time ti required for the velocity profiles to reach an 
approximately linear profile versus (a) <j> at fixed b„ = 0.0375 
and (b) inverse damping coefficient l/b n at fixed (f> = 0.90 for 
systems sheared at it = 0.364. In panel (b), the solid line has 
slope 0.14. 

has significant excursions above and below linearity at 
short times signals that the system is in the quasistatic 
flow regime. This behavior at small u is interesting but 
not the topic of this study. 

In Fig. we show that the time U required for the 
velocity profile to become linear increases with the ve- 
locity u of the shearing boundary. More precisely, ti — t s 
appears to scale as a power-law in u — ui 



(u - m) 







(6) 



for u > ui, where /3 rj 0.48±0.02. For u < uj, we find that 
the velocity profiles are nonlinear only for times t < t s , 
and are linear for all subsequent times. 

We also investigated the sensitivity of ti to changes in 
the packing fraction and damping coefficient. In Fig. [S] 
(a), we show that U increases sharply near random close 
packing rcp at small damping. This rise in ti is sup- 
pressed at large damping coefficients as shown in Fig. [3] 
The rapid rise in ti at least for undcrdamped systems at 
densities below random close packing again suggests that 
these systems are susceptible to the formation of strongly 
nonlinear velocity profiles. 

These results are surprising and raise an important 
question regarding the physical mechanism that is re- 
sponsible for the slow evolution of the velocity profiles. 
As a first step in addressing this question, we show be- 
low that granular temperature differences across the sys- 
tem give rise to nonlinear velocity profiles, and that, if a 
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FIG. 6: Ensemble averaged mean-square velocity fluctua- 
tions (Svy) in the shear gradient direction versus height y/L y 
measured from the stationary boundary for the same sys- 
tem shown in Fig. [5] The times shown are t = 520 (dashed 
line), 1000 (dotted line), 2000 (dot-dashed line), and 20000 
(solid). The mean-square velocity fluctuations become uni- 
form at long times. 

sufficiently large granular temperature difference can be 
maintained, nonlinear velocity profiles will be stable at 
long times. 



B. Combining vibration and shear: Stabilizing 
nonlinear velocity profiles at long times 

Recent experiments on sheared granular materials have 
shown that strongly nonlinear velocity profiles are accom- 
panied by spatially dependent granular temperature pro- 
files 0- Thus, an important question to ask is what role 
does the granular temperature play in determining the 
shape of the velocity profiles in sheared granular systems. 
Also, do these systems require a sufficiently large gran- 
ular temperature difference across the system to possess 
strongly nonlinear velocity profiles? Our results in Fig.[S] 
suggest that the shapes of the granular temperature and 
velocity profiles are strongly linked. In this figure, we 
show the time evolution of the mean-square velocity fluc- 
tuations in the shear-gradient direction, {Svy), following 
the initiation of shear. The sequence of times is identical 
to that shown in Fig. [21 At short times, there is a large 
difference in the mean-square velocity fluctuations be- 
tween the 'hot' shearing boundary and 'cold' stationary 
boundary, and the velocity profile is highly nonlinear. In 
contrast, at long times, the mean-square velocity fluctu- 
ations are uniform and the velocity profile is linear. Note 
that in the systems studied here, the granular tempera- 
ture difference across the system, AT, is roughly equal 
to the granular temperature near top boundary, since the 
mean-square velocity fluctuations are much smaller near 
the bottom stationary boundary. 

Figs. |3 and [S] show that large granular temperature 
differences and nonlinear velocity profiles occur together. 
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FIG. 7: The degree of nonlinearity (A(i)) y as a function 
of time for a system that is both sheared and vibrated with 
4> = 0.90, u = 0.0727, e = 0.92, and w = 0.1, 0.4, 0.8, 1.4, and 
1.8. The indexes 1 to 5 label the frequencies from smallest to 
largest. 



However, under steady shear, the granular temperature 
profiles become uniform and the velocity profiles become 
linear at long times. To further investigate the connec- 
tion between the granular temperature and velocity pro- 
files, we study systems that are both sheared and vi- 
brated, and are thus designed so that granular tempera- 
ture differences across the system can be maintained at 
long times. As discussed above in Sec. [H] in our second 
set of simulations we drive the top boundary at constant 
horizontal velocity u, while the location of each particle 
in the top or bottom boundary oscillates sinusoidally in 
time with fixed amplitude A and frequency u. The ver- 
tical vibrations cause the mean-square velocity fluctua- 
tions and packing fraction to become spatially nonuni- 
form with higher velocity fluctuations and lower packing 
fraction or dilatancy near the vibrated boundary. 

Fig. [7] clearly demonstrates that vertical vibration cou- 
pled with shear flow gives rise to stable nonlinear velocity 
profiles at long times. At low and also at high vibration 
frequencies, the degree of nonlinearity (A(t)) y of the ve- 
locity profile decays to small values at long times, as we 
found previously in Fig. [3] for the unvibrated systems. 
In contrast, at frequencies near lo c , {A(t)) y is nonzero 
at long times. We note that the longest times shown in 
Fig. [7| are at least 10 times longer than ti in the corre- 
sponding unvibrated system, confirming that these are 
indeed steady-state results. 

The results presented in Figs. |H1 (a) and (b) suggest 
that large and sustained granular temperature differences 
across the system and the resulting dilatancy are respon- 
sible for stable nonlinear velocity profiles. This figure 
shows that nonlinear velocity profiles occur when there 
are large granular temperature differences at lo = 0.4, 
0.8, and 1.4, while linear profiles are found when the ve- 
locity fluctuations are uniform at to = 0.1. Moreover, the 
degree of nonlinearity in the velocity profiles increases 
with the magnitude of the granular temperature differ- 
ence across the system. Fig. |S] (a) also shows the glass 
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FIG. 8: Time and ensemble averaged profiles for (a) velocity 
fluctuations {Svy) t , (b) packing fraction (4>)t, and (c) flow 
velocity (v x )t/u for a system at cf> = 0.90 and b n = 0.0375 
with the top boundary moving at u = 0.0727 and vibrated at 
ui = 0.1 (circles), 0.4 (squares), 0.8 (diamonds), 1.4 (upward 
triangles), and 1.8 (downward triangles). The dotted line in 
panel (a) shows the value of the glass transition temperature 
in a quiescent system at the same average packing fraction. 
The solid lines in panel (c) show numerical fits to the linear 
portions of the velocity profiles at lo = 0.4 and 1.4. 

transition temperature at the same average packing frac- 
tion; we discuss the relevance of this temperature in more 
detail in Sec.HVl|2^. 

The granular temperature difference across the system 
increases with lo for lo < lo* but decreases when lo > lo* 
with lo* k, lo c . The largest granular temperature differ- 
ence occurs near the natural frequency lo c of the repulsive 
spring interactions [H|. The decrease for lo > uo* occurs 
because particles adjacent to the vibrating boundary do 
not have enough time to react to the collision with the 
boundary before another collision occurs. As a result, vi- 
brations at large frequencies simply reduce the effective 
height of the system by the amplitude of the vibration 
but do not induce large granular temperature differences. 
We note that lo* does not appear to depend on the shear- 
ing velocity u. 

It is important to emphasize the fact that differences in 
the mean-square velocity fluctuations across the system, 
not the magnitude of the fluctuations themselves, are im- 
portant in determining the shape of the velocity profiles. 
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FIG. 9: Time and ensemble averaged profiles for (a) velocity 
fluctuations (5vy) t in the shear gradient direction, (b) packing 
fraction {(j))t, and (c) velocity {v x )t/u in the shear flow direc- 
tion for the same system in Fig. [8] except the bottom (not 
the top) boundary vibrates vertically at u = 0.1 (circles), 
0.4 (squares), 0.8 (diamonds), 1.4 (upward triangles), and 1.8 
(downward triangles). The solid lines in panel (c) show nu- 
merical fits to the linear portions of the velocity profiles at 
u = 0.4 and 1.4. 



For example, a system with the same interactions, av- 
erage density, and uniform temperature (Svy) > 10~ 3 
will possess a linear velocity profile when sheared over 
the same range of u. Similarly, the system vibrated at 
u> = 0.1 (shown as circles in Fig. |HJ) is in a glassy state 
with small but relatively uniform mean-square velocity 
fluctuations (Sv^) « 10~ 4 and also possesses a linear ve- 
locity profile. 

In systems that possess nonlinear velocity profiles, we 
find that the largest local shear rate does not occur 
equally likely at both boundaries. Instead, the portion of 
the system with the largest local shear rate always forms 
near the boundary with the largest (Svy) and resulting 
dilatancy. This is confirmed in Fig. [5] which shows re- 
sults in a system that is identical to that in Fig. |H1 except 
the bottom wall, not the top wall, is vibrated. The verti- 
cal vibrations induce a larger granular temperature near 
the bottom vibrated but unsheared boundary. Thus, the 
largest local shear rates occur near the bottom boundary, 
as shown for uj = 0.4, 0.8, and 1.4 in Fig. 0(c). 



IV. DISCUSSION 

In this section, we will focus on several aspects of our 
results in more detail. First, we find that nonlinear ve- 
locity profiles form only when the difference in the gran- 
ular temperature across the system exceeds a threshold 
value AT > AT$. Second, in systems with large gran- 
ular temperature differences AT > ATo at long times, 
the velocity profiles are linear near the 'cold' wall but 
highly nonlinear near the 'hot' wall. These profiles differ 
in shape from those found in the sheared but unvibrated 
systems . Finally, we point out that shear bands form 
when the average shear stress of the system (<J xy )t falls 
below the yield stress o~o required to initiate sustained 
flow in a static system. 

Figs. |H1 and [5] clearly show that differences in the mean- 
square velocity fluctuations across the system give rise to 
nonlinear velocity profiles. However, our results indicate 
that the granular temperature difference required to gen- 
erate a nonlinear velocity profile must exceed a threshold 
value. For example, systems with vibration frequency 
uj = 1.8 (downward triangles) in Figs. |H1 and El have rela- 
tively large granular temperature differences (slightly less 
than AT = 10~ 3 ), but possess nearly linear velocity pro- 
files. In contrast, systems with AT > 10~ 3 (for example 
uj = 0.4, 0.8, and 1.4) possess strongly nonlinear velocity 
profiles. 

The threshold granular temperature difference, which 
is roughly equal to the granular temperature near the 
vibrated boundary, appears to agree with the glass tran- 
sition temperature for a quiescent equilibrium system at 
the same average packing fraction [22]. This correspon- 
dence seems reasonable since in equilibrium systems the 
temperature must exceed the glass transition tempera- 
ture to cause large density fluctuations. Similarly, in 
sheared dissipative systems, it is difficult to create suffi- 
ciently large density gradients required for nonlinear ve- 
locity profiles if the granular temperature difference is 
below a threshold AT). 

We also find that the nonlinear velocity profiles that 
occur in the sheared and vibrated athermal systems have 
qualitatively different shapes compared to those found 
for the sheared but unvibrated systems. For example, 
the nonlinear velocity profiles for uj — 0.4 (squares) and 
1.4 (upward triangles) in Fig. [H] are composed of a lin- 
ear portion that extends from the bottom boundary at 
y = to y ~ 0.8, and a strongly nonlinear part near 
the top boundary. We note that in vibrated systems the 
crossover between the linear and highly sheared behav- 
ior occurs where the local packing fraction switches from 
uniform to nonuniform. 

Finally, we will address an interesting claim made in 
Ref. 0| that shear bands form in sheared glassy systems 
when the average shear stress in the system falls below 
the yield shear stress required to induce flow in a static 
state. In systems that form shear bands, shear flow is 
confined to a small portion of the system while the re- 
mainder of the system remains nearly static. In contrast, 
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FIG. 10: Time averaged shear stress ((T X y)t versus the bound- 
ary velocity u for a sheared but unvibrated system at (j> — 0.90 
and b n — 0.0375. The dotted line indicates the yield stress 
<To required to initiate sustained flow in a static system at the 
same packing fraction. 



all parts of the system flow when systems possess generic 
nonlinear velocity profiles. 

Does the constraint on the average shear stress guar- 
antee that shear bands will occur in the repulsive ather- 
mal systems studied here? In Fig. we compare the 
time-averaged shear stress {<J X y)t sheared at fixed u to 
the yield shear stress cr required to induce flow in an 
originally static unsheared state at the same <\>. We find 
that at small boundary velocities u, {<J xy )t is slightly be- 
low (To- However, even though the average shear stress is 
below the yield stress, the velocity profiles are linear for 
long times (as shown in Fig. |3J) and not highly localized 
or shear-banded. 

Thus, the condition {<7 xy )t < ctq alone does not en- 
sure that sheared repulsive athermal systems will form 
shear bands. An additional requirement must be satis- 
fied to stabilize shear bands at long times — the systems 
must possess sufficiently large granular temperature dif- 
ferences. It should be noted, however, that the differences 
between the average and the yield shear stress are small 
in the systems we studied and that shear stress fluctu- 
ations may inhibit the formation of shear bands. Thus, 
more work should be performed to verify the presented 
results. In fact, we are now attempting to determine the 
variables that set the difference between {o X y)t and ctq 
and whether this difference persists in the large system 
limit. 

Fig. ^2 shows that the average shear stress falls be- 
low the yield shear stress in the sheared and vibrated 
systems over a range of frequencies at u = 0.00727 and 
0.0727, but not at 0.727. We have also confirmed that the 
granular temperature differences in these systems satisfy 
AT > ATo over the range of frequencies 0.4 < ui < 1.4. 
Therefore, from the discussion above, one expects shear 
bands to form for the two higher shear rates, but not the 
lowest one. This prediction is confirmed in Fig. ^| At 
u = 0.727, nonlinear velocity profiles form, but shear is 
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FIG. 11: Time averaged shear stress {a~ xy )t for the sheared 
and vibrated system at <j> = 0.90 and b„ — 0.0375 as a function 
of u for u = 0.00727 (circles), 0.0727 (squares), and 0.727 
(upward triangles). The dotted line indicates the yield shear 
stress (7o for the static system at the same average density 
and dissipation. 

not highly localized into a shear band. In contrast, at 
u = 0.00727 there is a range of frequencies 0.4 < u> < 1.4 
over which shear bands form. 



V. CONCLUSIONS 

In this article we reported on recent simulations 
of model frictionless granular systems undergoing 
boundary-driven planar shear flow in 2D over a range 
of flow velocities and average densities, and for particles 
with varying degrees of inelasticity. These studies have 
produced several interesting and novel results that are 
relevant to a variety of jammed and glassy systems sub- 
jected to planar shear flow. First, we find that nonlinear 
velocity profiles are not stable at long times. Nonlinear 
velocity profiles form when the boundary velocity exceeds 
a characteristic speed set by the shear wave speed in the 
material, but they slowly evolve toward linear profiles at 
long times. In addition, the granular temperature and 
packing fraction profiles are initially spatially dependent 
but become uniform at long times. We measured the 
time ti required for the velocity profiles to become linear 
and for the granular temperature and density to become 
homogeneous throughout the system. We find that ti 
increases as u 5 at large u and increases strongly as <ft 
approaches rcp for nearly elastic systems. 

These results imply that sufficiently large and sus- 
tained granular temperature differences between the 'hot' 
and 'cold' boundaries are required to stabilize nonlinear 
velocity profiles at long times. We also studied systems in 
which vertical vibrations of the top or bottom boundary 
were superimposed onto planar shear flow to maintain 
granular temperature differences across the system. In 
the sheared and vibrated systems, we find that nonlin- 
ear velocity profiles are stable when the granular tern- 
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FIG. 12: Time averaged velocity profiles (v x (y)}t in the 
sheared and vibrated systems at <j> = 0.90 and fe n = 0.0375 for 
several vibration frequencies ui = 0.1 (circles), 0.4 (squares), 
0.8 (diamonds), 1.4 (upward triangles), and 1.8 (downward 
triangles) at (a) u = 0.727, (b) 0.0727, and (c) 0.00727. 



perature difference exceeds a threshold value ATo , which 
roughly corresponds to the glass transition temperature 
in an equilibrium system at the same average density. 
The nonlinear velocity profiles, however, differ in shape 
from those found previously for planar shear flow. The 
velocity profiles are composed of a linear part that ex- 
ists where the packing fraction is spatially uniform and a 
nonlinear portion that exists where the packing fraction 
varies strongly in space. Finally, we have shown that the 
nonlinear velocity profiles become highly localized when 
the average shear stress in the system is below the yield 
shear stress. 



VI. FUTURE DIRECTIONS 

Several additional studies are necessary to fully under- 
stand dense shear flows in granular systems, and these 
will be presented in future work |'25|. First, we intend 
to investigate the influence of dynamic and static fric- 
tion forces on the long-time stability of nonlinear velocity 
profiles. Preliminary results suggest that dynamic fric- 
tion is not sufficient to stabilize nonlinear velocity profiles 





FIG. 13: Snapshots from 2D systems undergoing Couette 
shear flow at (j> = 0.845 and b n = 0.0375. Panels (b), (c), and 
(d) differ from panel (a) by 2, 6, and 19 rotations, respectively. 
All particles with centers inside the inner ring at r = Ri move 
at fixed rotation rate Q, — 0.01 counter clockwise. All particles 
with centers outside the outer ring at r — R2 are stationary. 

in systems undergoing planar shear flow. In initial stud- 
ies with weak dissipation and dynamic friction | 26| near 
0rcp, we found that ti is similar to that for systems with 
dissipation and no dynamic friction. However, more ex- 
tensive studies of dynamic as well as static friction are 
required to make a definitive statement about the influ- 
ence of friction on the long-time stability of nonlinear 
velocity profiles 0, . 

Second, most experimental studies of velocity profiles 
in granular systems are performed in (angular) Couette, 
not planar shear cells. How does the geometry of the 
shear cell influence the velocity profiles? Are shear bands 
stable in frictionless, athermal systems undergoing Cou- 
ette shear flow? Fig. ^5] shows preliminary results from 
studies of 2D Couette shear flow in systems at 4> = 0.845, 
b n = 0.0375, and rotation rate — 0.01 in the counter 
clockwise direction. Panels (b), (c), and (d) show snap- 
shots of the system 2, 6, and 19 rotations after the initial 
configuration in panel (a). The snapshots provided in 
panels (a) - (c) reveal that the shear band in this system 
is approximately 5 — 8 small particle diameters wide since 
the two highlighted particles closest to the outer bound- 
ary do not rotate significantly even after 6 rotations of 
the inner boundary. A comparison of panels (c) and (d) 
shows that at long times the particles in the flowing re- 
gion are able to diffuse perpendicular to the boundaries. 
The mean-square velocity fluctuations (Si)g) in the tan- 
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FIG. 14: (a) Mean-square velocity fluctuations (5vg(r)) in the 
tangential direction and (b) tangential velocity ve(r) normal- 
ized by the speed ve(Ri) at the inner shearing boundary ver- 
sus the distance from the inner boundary (r — R\)/(R2 — Ri) 
in a 2D Couette cell with <j> = 0.845 and b n = 0.0375. The 
symbols indicate that time averages taken between the 20th 
and 80th (circles) or the 60th and 80th (triangles) rotations 
are nearly identical. 

gcntial direction and the tangential velocity vg normal- 
ized by the speed at the inner boundary are shown in 
Figs. El ( a ) an d (b) as a function of the distance from 
the inner boundary (r — i?i)/(i?2 — Ri)- Averages over 



varied numbers of rotations demonstrate that the non- 
linear velocity profiles are stable at long times. As we 
found in our studies of planar shear flow, nonlinear ve- 
locity profiles occur when the granular temperature is 
nonuniform. However, in Couette shear flows, the shear 
stress is also spatially dependent. Thus, the distinct con- 
tributions from nonuniform shear stress and nonuniform 
granular temperature need to be disentangled. In future 
work, we will determine whether nonlinear velocity pro- 
files persist when we vibrate the outer boundary to create 
a more uniform granular temperature profile. 

Finally, there have been several recent computational 
studies of effective temperatures defined from fluctuation 
dissipation relations, linear response theory, and elastic 
energy fluctuations in dense granular systems |27ll28l.l29| . 
These studies have shown that a consistent effective tem- 
perature can be defined for dense shear flows, i.e. the 
effective temperatures T c g obtained from the above defi- 
nitions agree with each other but T e g is much larger than 
the granular temperature of the system. Because T e g de- 
scribes fluctuations on long length and time scales, it is 
possible that the effective temperature will play a signif- 
icant role in determining the shape of velocity profiles. 
Thus, it is important to study the effective temperature 
as well as the granular temperature in sheared glassy and 
athermal systems. 
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